# SLE human patients
source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')

proj.nm <- 'mission/FPP/xiangya_sle_scRNA/'

late.kc <- c('Krt1','Krt10','Lor','Ivl','Tgm1','Flg') |>
  str_to_upper()

# mva pathway ----------
# ACAT >> HMGCS >> HMGCR >> MVK >> PMVK >> MVD >> IDI >> (FDPS,GGPS1)
kegg_mva <-
  c('Acat1','Acat2','Hmgcs1','Hmgcs2','Hmgcr','Mvk','Pmvk','Mvd','Idi1','Idi2','Fdps') |>
  str_to_upper()

mva_fct <- tibble(gene = kegg_mva, ordered = fct_inorder(kegg_mva))

# interested cytokines ------
key_cytokine <- c('Ccl20','Tslp','Flt3lg','Csf2','Tnf','Il6') |>
  str_to_upper()

# select only epidermis ----
sobj_epi <- read_rds('mission/FPP/xiangya_sle_scRNA/data/human_SLE_epi.rds')

sobj_epi |> dplyr::count(dataset)

sobj_epi <- sobj_epi |>
  quick_process_seurat(res = 1)

sobj_epi |>
  DimPlot(group.by = 'manual_main', cols = DiscretePalette(36),
          label = T, repel = T)

## annotate cell type ------
hpca <- celldex::HumanPrimaryCellAtlasData()

sobj_epi <- sobj_epi |>
  mark_cell_type_singler(hpca, new_label = 'hpca_main')

sobj_epi |>
  DotPlot(c('KRT1','KRT14','PMEL','MLANA'),cluster.idents = T)

sobj_epi <- sobj_epi |>
  mutate(manual_main = ifelse(hpca_main == 'Neurons', 'Melanocytes', hpca_main))

sobj_epi |> write_rds('mission/FPP/human_SLE_epi.rds')

## FIG: SLE epi cell type umap ==========
sobj_epi |>
  DimPlot(group.by = 'manual_main', cols = DiscretePalette(36),
          label = T, label.size = 2) +
  ggtitle('Cell types of epidermis') +
  theme_jpub +
  NoLegend()

publish_pdf('mission/FPP/figures/human_epi_celltype_umap.pdf')

## examine DC -------
### FIG: DC frac change SLE vs HC ==========
dc_conf <- sobj_epi |>
  dplyr::count(sle, DC = manual_main == 'DC') |>
  group_by(sle) |>
  calc_frac_conf_on_grouped_count() |>
  filter(DC)

sobj_epi@meta.data |>
  as_tibble() |>
  mutate(DC = ifelse(manual_main == 'DC', 'DC', 'nonDC')) |>
  test_on_grouped_count(sle, DC)

dc_conf |>
  ggplot(aes(sle, fraction, ymax = conf.high, ymin = conf.low, fill = sle)) +
  geom_col() +
  geom_errorbar(width = .5) +
  labs(x = 'Group', y = 'Fraction in all epidermis cells') +
  scale_fill_manual(values = c('skyblue','orange')) +
  theme_jpub

ggsave('figures/SLE-HC_epi_DC_frac.pdf',
       width = 50, height = 50, units = 'mm')

### FIG: alter DC frac ==========
kim_dc <- read_csv('mission/FPP/xiangya/kim_dc_count.csv')

kim5dc <- kim_dc |>
  filter(orig.ident == 'HC5kim') |>
  add_case(orig.ident = 'HC5kim', DC = 'non', fraction = .99278) |>
  mutate(DC = DC == 'DC')

dc_conf <- sobj_epi |>
  dplyr::count(orig.ident, DC = manual_main == 'DC') |>
  group_by(orig.ident) |>
  calc_frac_conf_on_grouped_count()

dc_conf |>
  bind_rows(kim5dc) |>
  arrange(orig.ident) |>
  mutate(group = str_extract(orig.ident, 'HC|SLE'),
         width = ifelse(group == 'HC', .3, .75),
         orig.ident = str_c(group, 2:15 %/% 2),
         cell_type = ifelse(DC, 'DC', 'other cells') |>
           fct_relevel('other cells')) |>
  ggplot(aes(orig.ident, fraction, fill = cell_type)) +
  geom_col() +
  facet_grid(~group, scales = 'free', space = 'free') +
  scale_fill_manual(values = c('grey','red')) +
  labs(x = 'Sample', y = 'Fraction in all cells', fill = 'Cell type') +
  theme_jpub +
  RotatedAxis()
  
publish_pdf('figures/SLE_new_DC_frac.pdf')

dc_conf |>
  bind_rows(kim5dc) |>
  arrange(orig.ident) |>
  mutate(group = str_extract(orig.ident, 'HC|SLE'),
         width = ifelse(group == 'HC', .3, .75),
         orig.ident = str_c(group, 2:15 %/% 2),
         cell_type = ifelse(DC, 'DC', 'other cells') |>
           fct_relevel('other cells')) |>
  filter(DC) |>
  write_csv('SLE_epidermis_DC_fraction.csv')

### FIG: DEG of SLE vs HC DC =========
dc_sle_deg <- sobj_epi |>
  filter(manual_main == 'DC') |>
  FindMarkers(group.by = 'sle',
              ident.1 = 'SLE',
              only.pos = T)

dc_sle_upgo <- dc_sle_deg |>
  filter(p_val_adj < .05 & avg_log2FC > 1) |>
  rownames() |>
  enrichGO(OrgDb = 'org.Hs.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           qvalueCutoff = .05,
           readable = T)

dc_sle_upgo |>
  enri_simplify() |>
  pluck('result') |>
  head(10) |>
  publish_enrichment()

publish_pdf('figures/SLE-HC_DC_upgo.pdf', width = 60)

kera_slevhc <- sobj_epi |>
  filter(manual_main == 'Keratinocytes') |>
  FindMarkers(group.by = 'sle',
              ident.1 = 'SLE') |>
  as_tibble(rownames = 'gene')

sobj_epi |>
  filter(manual_main == 'Keratinocytes') |>
  VlnPlot(features = key_cytokine, group.by = 'sle', pt.size = 0)

# only kera --------------------
sobj_keras <- sobj_epi |>
  filter(manual_main == 'Keratinocytes')

sobj_keras <- sobj_keras |>
  quick_process_seurat(res = .7, skip_norm = T)

sobj_keras <- sobj_keras |>
  filter(seurat_clusters != 15)

sobj_keras |> write_rds('mission/FPP/xiangya_sle_scRNA/data/human_SLE_kera.rds')

sobj_keras <- read_rds('mission/FPP/xiangya_sle_scRNA/data/human_SLE_kera.rds')

## FIG: kera clusters umap ============
sobj_keras |>
  filter(umap_1 < 10) |>
  DimPlot(cols = DiscretePalette(36), label = T) +
  theme_jpub +
  NoLegend()

publish_pdf('mission/FPP/figures/SLE_HC_kera_umap.pdf')

## FIG: kera clusters umap split by SLE/HC ==========
sobj_keras |>
  filter(umap_1 < 10) |>
  DimPlot(cols = DiscretePalette(36), label = T,
          split.by = 'sle', label.size = 2) +
  theme_jpub +
  NoLegend()

publish_pdf('mission/FPP/figures/SLE_HC_kera_umap_split.pdf')

### FIG: HC-TRPV3 SLE-CCL20 kera dotplot ===========
sobj_keras |>
  filter(sle != 'SLE') |>
  DotPlot('TRPV3') +
  theme_jpub +
  scale_color_gradient2(low = 'blue', high = 'red') +
  labs(x = 'Gene', y = 'Keratinocyte clusters')

g1 <- last_plot()

sobj_keras |>
  filter(sle == 'SLE') |>
  DotPlot(c('CCL20','TRPV3'), group.by = 'seurat_clusters') +
  theme_jpub +
  labs(x = 'Gene', y = 'Keratinocyte clusters',title = 'SLE group')

g2 <- last_plot()

g1 + NoLegend() + g2

publish_pdf('figures/SLE_kera_ccl20_bubble.pdf')

### FIG: HC TRPV3 & SLE CCL20 kera enrich dotplot ===========
v3_rich <- sobj_keras |>
  filter(sle != 'SLE') |>
  FindAllMarkers(features = 'TRPV3', logfc.threshold = 0, return.thresh = 1)

v3_rich |>
  ggplot(aes(gene, cluster, size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  labs(y = 'Keratinocyte clusters', x = 'Gene', title = 'HC') +
  scale_color_gradient2(low = 'blue', high = 'red') +
  scale_radius(range = c(.1,4)) +
  theme_jpub

publish_pdf('figures/HC-kera-TRPV3_enrich.pdf', width = 40)

er1 <- last_plot()

Idents(sobj_keras) <- 'seurat_clusters'

sle20_rich <- sobj_keras |>
  filter(sle == 'SLE') |>
  FindAllMarkers(features = 'CCL20', logfc.threshold = 0)

sle20_rich |>
  ggplot(aes(gene, cluster, size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  labs(y = '', x = 'Gene', title = 'SLE') +
  scale_color_gradient2(low = 'blue', high = 'red') +
  theme_jpub

er2 <- last_plot()

er1 + NoLegend() + er2

publish_pdf('figures/HC-kera-TRPV3_enrich.pdf', width = 40)

### FIG: key cytokine in SLE vs HC kera ====
sle_hc_kera_cytk <- sobj_keras |>
  FindMarkers(group.by = 'sle',
              ident.1 = 'SLE',
              features = key_cytokine,
              logfc.threshold = 0,
              min.pct = 0) |>
  as_tibble(rownames = 'gene')

sle_hc_kera_cytk |>
  mutate(p_val_adj = ifelse(p_val_adj == 0, 1e-200, p_val_adj)) |>
  ggplot(aes(y = gene, x = '', size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  labs(x = 'SLE-Keratinocyte', y = 'Gene') +
  scale_color_gradient2(low = 'blue', high = 'red') +
  scale_size(breaks = c(0,100,200)) +
  theme_jpub

ggsave('figures/SLE-HC_kera_cytk_bubble.pdf',
       width = 50, height = 50, units = 'mm')

## define TRPV3-hi kera ---------
sobj_keras <- sobj_keras |>
  mutate(trpv3h = ifelse(seurat_clusters %in% c(2,4), 'TRPV3-hi', 'TRPV3-lo'))

### FIG: trpv3h kera umap ===========
sobj_keras |>
  filter(umap_1 < 10) |>
  DimPlot(group.by = 'trpv3h') +
  labs(title = 'All keratinocytes') +
  theme_jpub

publish_pdf('figures/SLE_HC_TRPV3hi_kera_umap.pdf')

### FIG: trpv3h with late KC marker -------
sobj_keras |>
  DotPlot(c('TRPV3', late.kc, 'KRT5', 'KRT14'),
          cols = 'RdYlBu', cluster.idents = T) +
  labs(x = 'Gene', y = 'Keratinocyte clusters',
       title = 'TRPV3 and late KC markers in human skin KC')

sobj_keras |>
  DotPlot(c('TRPV3', late.kc, 'KRT5', 'KRT14'),
          cols = 'RdYlBu', cluster.idents = T, dot.scale = 2.5) +
  labs(x = 'Gene', y = 'Keratinocyte clusters',
       title = 'TRPV3 and late KC markers in human skin KC') +
  theme_jpub +
  RotatedAxis()

publish_pdf('mission/FPP/figures/KC_V3_late_marker_bubble.pdf', width = 60)

### SLE v3h vs v3l deg -------------
sle_v3hvl <- sobj_keras |>
  filter(sle == 'SLE') |>
  FindMarkers(group.by = 'trpv3h',
              ident.1 = 'TRPV3-hi',
              logfc.threshold = 0)

sle_v3hvl |>
  as_tibble(rownames = 'gene') |>
  write_csv('mission/FPP/hs_sle_v3hvl.csv')

### v3h SLE vs HC deg -------------
v3h_slevhc <- sobj_keras |>
  filter(trpv3h == 'TRPV3-hi') |>
  FindMarkers(group.by = 'sle',
              ident.1 = 'SLE')

v3h_slevhc |>
  as_tibble(rownames = 'gene') |>
  write_csv('mission/FPP/hs_v3h_sle-hc_deg.csv')

### v3l SLE vs HC deg -------------
v3l_slevhc <- sobj_keras |>
  filter(trpv3h == 'TRPV3-lo') |>
  FindMarkers(group.by = 'sle',
              ident.1 = 'SLE')

v3l_slevhc |>
  as_tibble(rownames = 'gene') |>
  write_csv('mission/FPP/xiangya_sle_scRNA/hs_v3l_sle-hc_deg.csv')

### HC v3h vs HC v3l -------
hc.v3hvl <- sobj_keras |>
  filter(sle != 'SLE') |>
  FindMarkers(group.by = 'trpv3h', ident.1 = 'TRPV3-hi') |>
  as_tibble(rownames = 'gene')
  
hc.v3hvl |>
  write_csv('mission/FPP/xiangya_sle_scRNA/hs_hc_v3hvl.deg.csv')

### FIG: v3h frac SLE vs HC =========
sobj_keras |>
  ggplot(aes(sle, fill = trpv3h)) +
  geom_bar(position = 'fill', width = .5) +
  labs(x = 'Group', y = 'Fraction in all keratinocytes', fill = 'Cell type') +
  theme_jpub

sobj_keras |>
  summarise(n(), .by = c(trpv3h,orig.ident)) |>
  pivot_wider(names_from = trpv3h, values_from = 'n()')

sobj_keras |>
  as_tibble() |>
  calc_frac_conf_on_grouped_count(sle, trpv3h) |>
  filter(trpv3h == 'TRPV3-hi') |>
  ggplot(aes(sle, fraction, ymin = conf.low, ymax = conf.high,
             color = sle)) +
  geom_col(fill = 'white') +
  geom_errorbar(width = .5) +
  theme_pubr() +
  scale_color_manual(values = c('blue','red')) +
  labs(x = 'Group', y = 'Fraction in all keratinocytes',
       title = 'Fraction of TRPV3-hi KC in skin',
       color = 'Group') +
  theme_jpub

publish_pdf('mission/FPP/figures/v3h_kc_frac.pdf')

sobj_keras |>
  as_tibble() |>
  test_on_grouped_count(sle, trpv3h)

## FIG: TRPV3 SLE-HC logfc violin ----------
trpv3h_meta <- sobj_keras |>
  as_tibble() |>
  select(.cell, sle, trpv3h)

trpv3h_meta <- sobj_keras |> 
  get_abundance_sc_long('TRPV3') |>
  left_join(trpv3h_meta)

v3pval <- tibble(group1 = 'SLE', group2 = 'HC',
                 y.position = 3.6, label = c('***','NS'),
                 trpv3h = c('TRPV3-hi','TRPV3-lo'))

trpv3h_meta |>
  ggplot(aes(sle, .abundance_RNA, fill = sle)) +
  geom_violin() +
  stat_pvalue_manual(data = v3pval, label = 'label',
                     label.size = 2, inherit.aes = F) +
  facet_wrap(~trpv3h, ncol = 1) +
  theme_pubr() +
  scale_fill_manual(values = c('blue','red')) +
  expand_limits(y = 4) +
  labs(x = 'Group', y = 'Normalized expression',
       title = 'TRPV3 expression', fill = 'Group')

trpv3h_meta |>
  ggplot(aes(sle, .abundance_RNA, fill = sle)) +
  geom_violin() +
  stat_pvalue_manual(data = v3pval, label = 'label',
                     label.size = 2, inherit.aes = F) +
  facet_wrap(~trpv3h, ncol = 1) +
  theme_pubr() +
  scale_fill_manual(values = c('blue','red')) +
  expand_limits(y = 4) +
  labs(x = 'Group', y = 'Normalized expression',
       title = 'TRPV3 expression', fill = 'Group') +
  theme_jpub

publish_pdf('figures/v3h_l_trpv3_logfc_sle_violin.pdf')

### FIG: all kera TRPV3 SLE-HC logfc =========
kera_14 <- sobj_keras$seurat_clusters |>
  unique()

kera_14_v3change <- kera_14 |>
map(\(x)sobj_keras |>
  FindMarkers(features = 'TRPV3',
              group.by = 'sle',
              ident.1 = 'SLE',
              subset.ident = x,
              logfc.threshold = 0),
  .progress = T) |>
  set_names(kera_14) |>
  list_rbind(names_to = 'cluster')

kera_14_v3change |>
  mutate(cluster = fct_inseq(cluster)) |>
  ggplot(aes(y = cluster, x = 'SLE vs HC TRPV3', size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  labs(x = '', y = 'Keratinocyte clusters') +
  scale_color_gradient2(low = 'blue', high = 'red') +
  theme_jpub

publish_pdf('figures/kera_trpv3_SLE-HC_logfc.pdf')

### FIG: v3h v3l TRPV3 SLE-HC logfc =========
v3_fct <- sobj_keras$trpv3h |>
  unique()

Idents(sobj_keras) <- 'trpv3h'

v3hl_v3change <- v3_fct |>
  map(\(x)sobj_keras |>
        FindMarkers(features = 'TRPV3',
                    group.by = 'sle',
                    ident.1 = 'SLE',
                    subset.ident = x,
                    logfc.threshold = 0),
      .progress = T) |>
  set_names(v3_fct) |>
  list_rbind(names_to = 'cluster')

v3hl_v3change |>
  ggplot(aes(y = cluster, x = 'SLE vs HC TRPV3', size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  labs(x = '', y = 'Keratinocyte clusters') +
  scale_color_gradient2(low = 'blue', high = 'red') +
  theme_jpub

publish_pdf('figures/kera_v3hl_trpv3_SLE-HC_logfc.pdf')

## SREBPs in V3h KC ---------
sobj_keras |>
  filter(sle != 'SLE') |>
  mutate(trpv3h = fct_relevel(trpv3h, 'TRPV3-lo')) |>
  bill.violin(c('SREBF1','SREBF2'), trpv3h) +
  labs(x = 'KC subtype', fill = 'KC subtype', y = 'Normalized RNA expression',
       title = 'SREBF1/2 expression in HC TRPV3-high/low keratinocytes')

sobj_keras |>
  filter(sle == 'SLE') |>
  mutate(trpv3h = fct_relevel(trpv3h, 'TRPV3-lo')) |>
  bill.violin(c('SREBF1','SREBF2'), trpv3h) +
  labs(x = 'KC subtype', fill = 'KC subtype', y = 'Normalized RNA expression',
       title = 'SREBF1/2 expression in SLE TRPV3-high/low keratinocytes')

sobj_keras |>
  filter(trpv3h == 'TRPV3-hi') |>
  #mutate(sle = fct_relevel(sle, 'HC')) |>
  bill.violin(c('SREBF1','SREBF2'), sle) +
  labs(x = 'Group', fill = 'Group', y = 'Normalized RNA expression',
       title = 'SREBF1/2 expression in SLE vs HC TRPV3-high keratinocytes')

## TRP family -----------
trp_gene <- rownames(sobj_keras) |> str_subset('TRP[CMV].$') |> str_sort()

sobj_keras |> filter(sle == 'HC') |>
  DotPlot(trp_gene, group.by = 'manual_main', cols = c('lightgrey','red')) +
  RotatedAxis()

last_plot() |>
  pluck('data') |>
  mutate(gene = fct_reorder(features.plot, avg.exp)) |>
  ggplot(aes(avg.exp, gene, fill = pct.exp)) +
  geom_col() +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  scale_fill_distiller(palette = 'Reds', direction = 1) +
  labs(x = 'Average expression', fill = 'Percent expressed',
       title = 'TRP family gene expression in HC human KC') +
  theme_jpub

publish_source_plot('hc.trp.family.expr.kc', width = 60)

## UPR chaperone ------
erupr.chap <-
  read_csv('mission/FPP/UPR.associated.chaperone.gene.csv')

sobj_keras |>
  filter(sle == 'SLE') |>
  DotPlot(erupr.chap$SYMBOL, group.by = 'trpv3h', scale = F,
          cols = c('white','red'))

last_plot() |>
  pluck('data') |>
  as_tibble() |>
  mutate(features.plot = fct_reorder(features.plot, avg.exp)) |>
  arrange(features.plot) |>
  write_csv('mission/FPP/pub_source_data/sle.v3kc.chaperone.expr.csv')

# back to all epi ---------
v3h_cell <- sobj_keras |>
  filter(trpv3h == 'TRPV3-hi') |>
  pull(.cell)

v3l_cell <- sobj_keras |>
  filter(trpv3h == 'TRPV3-lo') |>
  pull(.cell)

sobj_epi <- sobj_epi |>
  mutate(bill_fine = case_when(.cell %in% v3h_cell ~ 'TRPV3-hi-KC',
                               .cell %in% v3l_cell ~ 'TRPV3-lo-KC',
                               str_detect(manual_main, 'Kera') ~ 'badkera',
                               .default = manual_main))

sobj_epi <- sobj_epi |>
  mutate(bill_fine = case_match(bill_fine, 'TRPV3-hi-Keratinocytes' ~ 'TRPV3-hi-KC',
                                'TRPV3-lo-Keratinocytes' ~ 'TRPV3-lo-KC',
                               .default = manual_main))

sobj_epi <- sobj_epi |>
  filter(bill_fine != 'badkera')

Idents(sobj_epi) <- 'bill_fine'

sobj_epi |> write_rds('mission/FPP/xiangya_sle_scRNA/data/human_SLE_epi.rds')

sobj_epi <- read_rds('mission/FPP/xiangya_sle_scRNA/data/human_SLE_epi.rds')

## MVA --------------
### FIG: mva expr in SLE or HC epi =========
sobj_epi |>
  filter(sle == 'SLE') |>
  DotPlot(kegg_mva, group.by = 'bill_fine') +
  labs(x = 'Gene', y = 'Cell type',
       title = 'MVA pathway in SLE human skin') +
  scale_color_gradient2(low = 'blue', high = 'red') +
  scale_size(range = c(.1,8)) +
  RotatedAxis()

sobj_epi |>
  filter(sle != 'SLE') |>
  DotPlot(kegg_mva, group.by = 'manual_main') +
  labs(x = 'Gene', y = 'Cell type',
       title = 'MVA pathway in HC human skin') +
  scale_size(range = c(.1,8)) +
  RotatedAxis()

hc_epi_mva <- last_plot() |>
  pluck('data')

hc_epi_mva |>
  mutate(id = fct_relevel(id, 'Keratinocytes', after = Inf)) |>
  BubblePlot() +
  theme_bw(base_size = 6, base_family = 'ArialMT') +
  RotatedAxis() +
  theme_jpub

publish_source_plot('HC.skin.mva.expr', width = 70)

sobj_epi |>
  filter(sle == 'SLE') |>
  DotPlot(kegg_mva, group.by = 'bill_fine', dot.scale = 3,
          cols = 'RdYlBu') +
  labs(x = 'Gene', y = 'Cell type',
       title = 'MVA pathway in SLE human skin') +
  theme_jpub +
  RotatedAxis()

publish_pdf('mission/FPP/figures/SLE_epiderm_mva_bubble.pdf', width = 70)

sobj_epi |>
  filter(sle != 'SLE') |>
  DotPlot(kegg_mva, group.by = 'bill_fine', dot.scale = 2,
          cols = 'RdYlBu') +
  labs(x = 'Gene', y = 'Cell type',
       title = 'MVA pathway in HC human skin') +
  theme_jpub +
  RotatedAxis()

publish_pdf('mission/FPP/figures/HC_epi_mva_bubble.pdf', width = 70)

last_plot() |>
  pluck('data') |>
  write_csv('mission/FPP/pub_source_data/HC.human.skin.mva.expr.csv')

mva_fct$ordered

### FIG: mva enrich in SLE epi ==========
epi_mva_rich <- sobj_epi |>
  filter(sle == 'SLE') |>
  FindAllMarkers(features = kegg_mva,
                 logfc.threshold = 0,min.pct = 0) |>
  left_join(mva_fct) |>
  mutate(p_val_adj = ifelse(p_val_adj == 0, 1e-300, p_val_adj))

epi_mva_rich |>
  ggplot(aes(ordered, as.character(cluster), size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  labs(x = 'Gene', y = 'Cell type') +
  scale_color_gradient2(low = 'blue', high = 'red') +
  scale_size(range = c(.1,9)) +
  theme_pubr(legend = 'right') +
  RotatedAxis()

epi_mva_rich |>
  ggplot(aes(ordered, as.character(cluster), size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  labs(x = 'Gene', y = 'Cell type') +
  scale_color_gradient2(low = 'blue', high = 'red') +
  scale_radius(range = c(.1,3)) +
  theme_jpub +
  RotatedAxis()

publish_pdf('figures/SLE_epiderm_mva_enrich.pdf', width = 70)

### FIG: SLE v3h vs v3l mva logfc over other epi dotplot ==========
epi_fine <- unique(sobj_epi$bill_fine) |>
  str_subset('E|B', negate = T)

epi_fine_fcs <- sobj_epi |>
  filter(!str_detect(bill_fine, 'E|B')) |>
  FindMarkersAcrossVar(split.by = 'bill_fine', group.by = 'sle',
                       ident.1 = 'SLE', logfc.threshold = 0)

epi_fine_fcs |>
  write_csv('mission/FPP/xiangya_sle_scRNA/epi.slevhc.all.logfc.csv')

epi_fine_fcs <- epi_fine |>
  map(\(x)FindMarkers(sobj_epi,
                      group.by = 'sle', ident.1 = 'SLE', subset.ident = x,
                      features = kegg_mva,
                      logfc.threshold = 0, min.pct = 0) |>
        as_tibble(rownames = 'gene'),
      .progress = T) |>
  setNames(epi_fine) |>
  list_rbind(names_to = 'celltype')

epi_fine_fcs |>
  left_join(mva_fct) |>
  ggplot(aes(x = ordered, y = celltype, size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  labs(y = 'Cell type', x = 'Gene',
       title = 'Differential expression of MVA pathway in SLE vs HC skin') +
  scale_color_gradient2(low = 'blue', high = 'red') +
  theme_pubr(legend = 'right', x.text.angle = 45)

epi_fine_fcs |>
  left_join(mva_fct) |>
  mutate(avg_log2FC = ifelse(avg_log2FC < -2 & p_val_adj > .05, -1.72, avg_log2FC)) |>
  ggplot(aes(x = ordered, y = celltype, size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  labs(y = 'Cell type', x = 'Gene',
       title = 'Differential expression of MVA pathway in SLE vs HC skin') +
  scale_color_distiller(palette = 'RdYlBu') +
  scale_size(breaks = c(0,100,200), range = c(0,4)) +
  theme_jpub +
  RotatedAxis()

publish_pdf('mission/FPP/figures/epiderm_SLE-HC_mva_logfc.pdf', width = 70)

### mva module score -----------
sobj_epi <- sobj_epi |>
  AddModuleScore(features = list(kegg_mva), name = 'MVA_module')

sobj_epi |>
  VlnPlot('MVA_module1', group.by = 'bill_fine', pt.size = 0) +
  stat_summary(geom = 'crossbar', fun = mean)

mva_score <- sobj_epi |>
  select(bill_fine, MVA_module1, sle) |>
  mutate(normalized_score = (MVA_module1 - mean(MVA_module1))/sd(MVA_module1))

### FIG: normalized mva score in cell types -----------
mva_score |>
  ggplot(aes(bill_fine, normalized_score, fill = bill_fine)) +
  stat_mean(geom = 'col') +
  stat_summary(geom = 'errorbar', fun.data = mean_cl_normal, width = .3) +
  facet_wrap(~sle) +
  labs(x = 'Cell type', y = 'Normalized module score') +
  theme_jpub +
  NoLegend() +
  RotatedAxis()

publish_pdf('mission/FPP/figures/epiderm_SLE-HC_mva_module.pdf', width = 70)

mva_score |>
  summarise(mean_score = mean(normalized_score), .by = c(sle, bill_fine)) |>
  ggplot(aes(sle, mean_score, color = bill_fine, group = bill_fine)) +
  geom_point() +
  geom_line()

mva_score |>
  pivot_wider(values_from = MVA_module1, names_from = sle,
              id_cols = -normalized_score, values_fn = mean)
  
### gsea of mva --------------
sle_epitype_marker <- sobj_epi |>
  filter(sle == 'SLE') |>
  FindAllMarkers(group.by = 'bill_fine') |>
  filter(p_val_adj < .05) |>
  as_tibble()

proj.nm <- 'mission/FPP/xiangya_sle_scRNA/'

sle_epitype_marker |>
  write_source_csv('SLE.epi.allmarker')

sle_epitype_marker$cluster |>
  unique() |>
  walk(\(x)sle_epitype_marker |>
         filter(p_val_adj < .05, avg_log2FC > 1, cluster == x) |>
         select(gene) |>
         write_source_csv(str_c('SLE_allmarker_', x), col_names = F))

hc_epitype_marker <- sobj_epi |>
  filter(sle != 'SLE') |>
  FindAllMarkers(group.by = 'bill_fine') |>
  filter(p_val_adj < .05) |>
  as_tibble()

hc_epitype_marker |>
  write_source_csv('HC.epi.allmarker')

hc_epitype_marker$cluster |>
  unique() |>
  walk(\(x)hc_epitype_marker |>
         filter(p_val_adj < .05, avg_log2FC > 1, cluster == x) |>
         select(gene) |>
         write_source_csv(str_c('HC_allmarker_', x), col_names = F))

kc_epitype_marker <- sobj_epi |>
  FindMarkersAcrossVar(split.by = 'sle', group.by = 'manual_main',
                       ident.1 = 'Keratinocytes')

kc_epitype_marker$cluster |>
  unique() |>
  walk(\(x)kc_epitype_marker |>
         filter(p_val_adj < .05, avg_log2FC > 1, cluster == x) |>
         select(gene) |>
         write_source_csv(str_c(x, '_allmarker_Keratinocytes'), col_names = F))

sle_epitype_marker |>
  filter(str_detect(cluster, 'TRPV3-hi')) |>
  select(gene, avg_log2FC) |>
  arrange(avg_log2FC) |>
  write_tsv('v3h_epi.rnk')

mva_term2gene <- tibble(term = 'MVA metabolism', gene = kegg_mva)

term2name <- tibble(term = 'MVA metabolism',name  = 'MVA metabolism')

mva_gsea <- v3h_deg$avg_log2FC |>
  set_names(v3h_deg$gene) |>
  sort(decreasing = T) |>
  clusterProfiler::GSEA(TERM2GENE = mva_term2gene, pvalueCutoff = 1,pAdjustMethod = 'none')

mva_gsea |> enrichplot::gseaplot2('MVA metabolism', pvalue_table = T)

### GSVA of MVA --------------
gsva_mva <- list('MVA metabolism' = kegg_mva)

system.time(gsva_mva_res <- sobj_epi |>
  filter(str_detect(bill_fine, 'TRPV3-hi')) |>
  GetAssayData() |>
  as.matrix() |>
  gsvaParam(gsva_mva, minSize = 3) |>
  gsva(BPPARAM = SnowParam(progressbar = T)))

v3h_sle_meta <- sobj_epi |>
  filter(str_detect(bill_fine, 'TRPV3-hi')) |>
  pull(sle)

modelmtx <- model.matrix(~ v3h_sle_meta)

modelmtx |> head()

library(limma)

gsva_mva_res |>
  lmFit(modelmtx) |>
  eBayes() |>
  topTable(coef = 2)

## cytokine -------------
### FIG: cytokine expr in HC/SLE epi =========
sobj_epi |>
  filter(sle != 'SLE') |>
  DotPlot(key_cytokine, group.by = 'bill_fine') +
  labs(x = 'Gene', y = 'Cell type',
       title = 'Antibody-promoting factors in HC human skin') +
  scale_color_gradient2(low = 'blue', high = 'red') +
  RotatedAxis()

sobj_epi |>
  filter(sle != 'SLE') |>
  DotPlot(key_cytokine, group.by = 'bill_fine', dot.scale = 3,
          cols = 'RdYlBu') +
  labs(x = 'Gene', y = 'Cell type',
       title = 'Antibody-promoting factors in HC human skin') +
  theme_jpub +
  RotatedAxis()

publish_pdf('mission/FPP/figures/HC_epi_cytk_bubble.pdf', width = 60)

sobj_epi |>
  filter(sle == 'SLE') |>
  DotPlot(key_cytokine, group.by = 'bill_fine') +
  labs(x = 'Gene', y = 'Cell type',
       title = 'Antibody-promoting factors in SLE human skin') +
  scale_color_gradient2(low = 'blue', high = 'red') +
  RotatedAxis()

sobj_epi |>
  filter(sle == 'SLE') |>
  DotPlot(key_cytokine, group.by = 'bill_fine', dot.scale = 3,
          cols = 'RdYlBu') +
  labs(x = 'Gene', y = 'Cell type',
       title = 'Antibody-promoting factors in SLE human skin') +
  theme_jpub +
  RotatedAxis()

publish_pdf('mission/FPP/figures/SLE_epiderm_cytk_bubble.pdf', width = 60)

### CCL20 tpm percent in cell types ---------
ccl20.tpm.hep <- sobj_epi |>
  get_abundance_sc_wide('CCL20') |>
  mutate(CCL20 = expm1(CCL20)) |>
  left_join(x = sobj_epi, y = _) |>
  summarise(sum = sum(CCL20), .by = c(bill_fine, sle))

ccl20.tpm.hep |>
  filter(sle != 'HC') |>
  mutate(percent = sum / sum(sum))

### FIG: cytokine enrich in all SLE epi ============
Idents(sobj_epi) <- 'bill_fine'

epi_cytk_rich <- sobj_epi |>
  filter(sle == 'SLE') |>
  FindAllMarkers(features = key_cytokine,
              logfc.threshold = 0)

epi_cytk_rich |>
  mutate(p_val_adj = ifelse(p_val_adj == 0, 1e-300, p_val_adj)) |>
  ggplot(aes(gene, as.character(cluster), size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  labs(x = 'Gene', y = 'Cell type') +
  scale_color_gradient2(low = 'blue', high = 'red') +
  theme_jpub +
  scale_size(range = c(.1,2.5)) +
  RotatedAxis()

publish_pdf('figures/SLE_epiderm_cytk_enrich.pdf', width = 60)

### HC-SLE cytokine logfc in epi ==========
epi_fine <- unique(sobj_epi$bill_fine) |>
  str_subset('E|B', negate = T)

epi_fine_fcs <- epi_fine |>
  map(\(x)FindMarkers(sobj_epi,
                      group.by = 'sle', ident.1 = 'SLE', subset.ident = x,
                      features = key_cytokine,
                      logfc.threshold = 0, min.pct = 0) |>
        as_tibble(rownames = 'gene'),
      .progress = T) |>
  setNames(epi_fine) |>
  list_rbind(names_to = 'celltype')

sig_limit <- epi_fine_fcs |>
  filter(p_val_adj < .01) |>
  pull(avg_log2FC) |>
  max()

epi_fine_fcs %<>%
  mutate(avg_log2FC = ifelse(p_val_adj > .05 & abs(avg_log2FC) > sig_limit,
                             avg_log2FC / abs(avg_log2FC) * sig_limit, avg_log2FC))

epi_fine_fcs |>
  ggplot(aes(x = gene, y = celltype, size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  labs(y = 'Cell type', x = 'Gene',
       title = 'Differential expression of cytokine in SLE vs HC skin') +
  scale_color_gradient2(low = 'blue', high = 'red') +
  theme_pubr(legend = 'right')

epi_fine_fcs |>
  ggplot(aes(x = gene, y = celltype, size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  labs(y = 'Cell type', x = 'Gene',
       title = 'Differential expression of cytokine in SLE vs HC skin') +
  scale_color_distiller(palette = 'RdYlBu') +
  scale_size(breaks = c(0,100,200), range = c(.1,3)) +
  theme_jpub +
  RotatedAxis()

publish_pdf('mission/FPP/figures/epiderm_SLE-HC_cytk_logfc.pdf', width = 60)

# CXCL13 for plasma cell ----------
sobj_epi |> DotPlot2d('CXCL13', bill_fine, sle) +
  RotatedAxis() +
  labs(title = 'CXCL13 expression in SLE human epidermis',
       x = 'Cell type', y = 'Group') +
  scale_color_distiller(palette = 'Reds', direction = 1)

g1 <- last_plot()

slevhc.dc.t.fc <- sobj_epi |>
  filter(str_detect(bill_fine, 'DC|T|Mela')) |>
  FindMarkersAcrossVar(split.by = 'bill_fine', group.by = 'sle',
                       ident.1 = 'SLE')

slevhc.dc.t.fc |>
  filter(gene == 'CXCL13') |>
  ggplot(aes(cluster, y = '', color = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point() +
  theme_pubr(legend = 'right', x.text.angle = 45) +
  scale_color_gradient2(low = 'blue', high = 'red') +
  labs(title = 'CXCL13 fold change: SLE vs HC',
       x = 'Cell type', y = '')

g2 <- last_plot()

g1 / g2

sobj_epi |>
  filter(bill_fine == 'DC') |>
  FindMarkers(group.by = 'sle', features = 'CXCL13',
              ident.1 = 'SLE')

# CCR6+ T cell --------
sobj_epi$bill_fine |> unique()

sle.all.frac <- sobj_epi |>
  calc_frac_conf_on_grouped_count(sle, bill_fine)

sle.all.frac |>
  ggplot(aes(sle, fraction, fill = sle)) +
  geom_col() +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = .3) +
  facet_wrap(~bill_fine, scales = 'free_y')

sobj.t <- sobj_epi |>
  filter(bill_fine == 'T_cells')

sobj.t |>
  FindMarkers(features = 'CCR6', group.by = 'sle', ident.1 = 'SLE')

sobj.t |>
  bill.violin(features = 'CCR6', group.by = sle)

sobj.t |> DimPlot(group.by = 'seurat_clusters', split.by = 'sle')

sobj.t |> DotPlot(group.by = 'seurat_clusters', features = 'CCR6')

sobj.t %<>% quick_process_seurat(skip_norm = T)

monaco <- celldex::MonacoImmuneData()

sobj.t %<>% mark_cell_type_singler(monaco, fine_label = T,
                                   new_label = 'monaco.fine')

sobj.t |> DimPlot(group.by = 'monaco.fine')

sobj.t |>
  calc_frac_conf_on_grouped_count(sle, monaco.fine) |>
  ggplot(aes(sle, fraction)) +
  geom_col() +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = .3) +
  facet_wrap(~monaco.fine, scales = 'free_y')

sobj.t |>
  FindAllMarkers(features = 'CCR6', only.pos = T)

sobj.t |>
  DotPlot(seurat_markers)

sobj.t |>
  DotPlot('CCR6',group.by = 'sle', split.by = 'monaco.fine',
          cols = DiscretePalette(4))

sobj.t |>
  mutate(subgroup = str_c(seurat_clusters, '_', sle)) |>
  DotPlot('CCR6', group.by = 'subgroup')

# ATP receptor ------
p2r.listh <- rownames(sobj_keras) |> str_subset('^P2')

p2r.listh |> sort()

sobj_epi |>
  filter(sle == 'HC') |>
  DotPlot(p2r.listh, cluster.idents = T) +
  RotatedAxis() +
  ggtitle('ATP receptor expression in HC human skin')

# TLRs --------
tlr.listh <- rownames(sobj_epi) |>
  str_subset('^TLR\\d+$') |>
  str_sort(numeric = T)

sobj_epi |>
  filter(sle == 'HC') |>
  DotPlot(tlr.listh, cols = 'RdYlBu', dot.scale = 3) +
  labs(title = 'TLRs in HC human skin',
       x = 'Gene', y = 'Cell type') +
  theme_jpub +
  RotatedAxis()

publish_pdf('mission/FPP/figures/tlr.in.hc.kc.pdf', width = 60, height = 40)

g1 <- last_plot()

g1$data |>
  write_csv('mission/FPP/pub_source_data/fig.A.TLR.in.HC.human.skin.csv')

# IFN receptor? -------
ifnr.listh <- rownames(sobj_epi) |> str_subset('^IFN.R')

sobj_epi |>
  filter(sle == 'HC') |>
  DotPlot(c(ifnr.listh,'IL10RB')) +
  ggtitle('IFN receptor expression in HC human skin')

ifnl.listh <- rownames(sobj_epi) |> str_subset('IFNL\\d')

sobj_epi |>
  filter(sle != 'HC') |>
  DotPlot(ifnl.listh) +
  ggtitle('IFN lamdba expression in SLE human skin')

sobj_epi |>
  mutate(subgroup = str_c(bill_fine, '_', sle)) |>
  DotPlot(infl.listh, group.by = 'subgroup')

sobj_epi |>
  VlnPlot('IFNL1', group.by = 'sle')

# NLRP3 inflammasome pathway -----------
nlrp3i.listh <- c('NLRP1','NLRP3','PYCARD','CASP1','IL1B','IL18','IL33',
                  'NLRP6','NLRP7','NLRP12','NLRC4','AIM2','MEFV','IFI16')

sobj_epi |>
  filter(bill_fine != 'B_cell') |>
  mutate(subgroup = str_c(bill_fine, '-', sle)) |>
  DotPlot(nlrp3i.listh, group.by = 'subgroup') +
  ggtitle('Inflammasome pathway in SLE skin') +
  RotatedAxis()

# FDPS in hpca --------
fdps.hpca <- hpca |> filter(.feature == 'FDPS')
fdps.hpca |>
  mutate(cell.type = fct_reorder(label.main, logcounts),
         type = ifelse(str_detect(label.main, 'HSC|Kera|tem'), 'good','bad')) |>
  ggplot(aes(y = cell.type, x = logcounts, color = type)) +
  geom_boxplot() +
  labs_pubr() +
  scale_color_manual(values = c('black','red')) +
  labs(title = 'FDPS RNA expression in Human Primary Cell Altas') +
  NoLegend()

# ER stress --------
hc.v3hvl.upkegg <- read_csv('mission/FPP/xiangya_sle_scRNA/hc.v3hvl.upkegg.csv')
hc.v3hvl.upkegg |> DT::datatable()

hc.v3hvl.upgo <- read_csv('mission/FPP/xiangya_sle_scRNA/hc.v3hvl.upgo.csv')
hc.v3hvl.upgo |> DT::datatable()

kc.slevhc.upgo <- read_csv('mission/FPP/xiangya_sle_scRNA/kc.slevhc.upgo.csv')
kc.slevhc.upgo |> DT::datatable()

# ROS metabolism ----------
kc.slevhc.up.ros.hit9 <-
kc.slevhc.upgo |>
  filter(ID == 'GO:0072593') |>
  separate_longer_delim(geneID, '/') |>
  mutate(hit = anno_pmc_hits(str_glue('{geneID} reactive oxygen'))) |>
  select(geneID, hit) |>
  slice_max(hit, n = 9)

sobj_epi |>
  filter(sle == 'SLE') |>
  DotPlot(kc.slevhc.up.ros.hit9$geneID, group.by = 'bill_fine',
          cols = 'RdYlBu') +
  RotatedAxis()

sobj_epi |>
  filter(sle != 'SLE') |>
  DotPlot(kc.slevhc.up.ros.hit9$geneID, group.by = 'bill_fine',
          cols = 'RdYlBu') +
  RotatedAxis()

ers.slevhc <- epi_fine |>
  map(\(x)sobj_epi |> FindMarkers(group.by = 'sle', ident.1 = 'SLE',
                              subset.ident = x,
                              features = v3h.upers.hit9$gene.hs) |>
        as_tibble(rownames = 'gene'),
      .progress = T) |>
  setNames(epi_fine) |>
  list_rbind(names_to = 'cluster')

ers.slevhc |>
  ggplot(aes(x = gene, y = cluster, size = -log10(p_val_adj), color = avg_log2FC)) +
  geom_point() +
  labs(y = 'Cell type', x = 'Gene',
       title = 'Differential expression of ER stress pathway in SLE vs HC skin') +
  scale_color_distiller(palette = 'RdYlBu') +
  theme_pubr(legend = 'right')

## ADRB2 ADRA1A ---------
sobj_keras |>
  filter(sle == 'HC') |>
  DotPlot(c('ADRB2','ADRA1A'))

adrs <- last_plot() |>
  pluck('data')

adrs |>
  pivot_wider(id_cols = id, names_from = features.plot, values_from = avg.exp) |>
  ggplot(aes(ADRB2, ADRA1A)) +
  geom_vline(xintercept = 0) +
  geom_hline(yintercept = 0) +
  geom_point() +
  labs(title = 'Human healthy KC') +
  theme_bw()

